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A Thermodynamic  Surface  for  Water: 

The  Formulation  and  Computer  Programs 

by 

Lester  Haar,  John  s.  Gallagher  and  George  S.  Kell 
1.  Introduction 

This  report  has  been  prepared  to  provide  computational  details  and 
computer  programs  for  the  evaluation  of  the  thermodynamic  properties  of 
water  and  steam.  It  is  based  on  an  extensive  correlation  of  thermodynamic 
properties.1  Included  herein  are  the  basic  equations  for  the  description 
of  the  thermodynamic  properties  of  water  and  detailed  listings  of  computer 
programs  based  on  these  equations.  With  this  information,  thermodynamic 
properties  can  be  calculated  from  the  freezing  point  to  '2500  K in 
temperature  and  from  the  dilute  gas  to  well  in  excess  of  1 GPa  (10,000 
bars)  in  pressure  for  liquid  and  gaseous  states  for  undissociated  water. 

In  reducing  the  data  the  molecular  weight  of  ordinary  water  has 
been  taken  as  18.0152  g mol"1  ,2  the  universal  gas  constant  as  8.31441 
J mol-1  K"1,  and  hence  the  specific  gas  constant  as  0.461522  J g"1  K"1.3 
Absolute  temperatures  on  the  International  Practical  Temperature  Scale 
of  1968  have  been  used  as  the  realization  of  the  thermodynamic  scale. 

This  report  contains  five  sections  and  2 appendices:  A brief 
discussion  of  the  structure  of  the  thermodynamic  surface  for  water  is 
given  in  section  2,  and  in  section  3 are  the  equations,  parameters  and 
constants  that  define  the  surface.  Section  4 contains  the  thermodynamic 
relations  used  in  the  computer  program.  It  also  contains  a discussion 
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of  the  general  organization  and  structure  of  the  computer  program. 

Section  5 contains  statements  regarding  the  accuracy  and  range  of 
validity  of  our  surface.  In  Appendix  A are  the  computer  listings  and  a 
quide  to  their  use,  and,  lastly,  in  appendix  B,  is  a printout  of  calculations 
obtained  with  the  program,  so  that  the  user  can  check  his  results. 

2.  Structure  of  the  thermodynamic  surface 

The  thermodynamic  surface  is  composed  of  three  parts:  (1)  The 
first,  referred  to  as  the  base  function,  is  obtained  from  a theoretical 
equation  of  state.4-12  It  properly  describes  the  low-density  vapor,  the 
high  temperature  behavior  at  all  densities,  and  the  dense  fluid  at  all 
temperatures.  Except  for  a large  region  around  the  critical  point,  the 
base  function  yields  results  that  are  at  least  in  semi -quantitative 
accord  with  data.  (2)  The  second,  referred  to  as  the  residual  function, 
yields  corrections  to  the  base  function.  These  corrections  are  small  in 
regions  where  the  base  function  is  in  good  accord  with  data,  and  in 
regions  beyond  the  range  of  the  data.  The  contributions  to  the  pressure 
from  the  base  function  and  the  residual  function  are  readily  integrable 
in  closed  form  to  yield  their  respective  contributions  to  the  Helmholtz 
function.  (3)  The  third  includes  the  thermodynamic  properties  for  the 
ideal  gas  state.  These  have  been  reported  recently  by  Woolley13  as 
part  of  a detailed  analysis  of  the  rotation-vibration  structure  of  the 
water  molecule. 

A major  part  of  the  work  was  given  to  the  derivation  of  the  40 
terms  that  form  the  residual  function.  The  first  36  terms  were  used  in 
a global,  least  squares  fit  to  data.  Each  of  these  terms  yields  important 
contributions  over  wide  ranges  of  the  independent  variables.  Following 
this,  small  improvements  were  made  by  adding  3 terms  that  contribute 
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only  in  the  immediate  neighborhood  of  the  critical  point,  and  a single 
term  that  contributes  only  in  the  region  of  high  pressures  and  low 
temperatures.  Except  in  these  very  limited  regions,  the  residual 
function  is  given  by  the  first  36  terms.  A discussion  of  the  thermodynamic 
surface  obtained  with  the  residual  function  so  restricted  (terms  1-36) 
is  given  in  reference  [12]. 

The  three  parts  are  combined  to  obtain  the  expression  for  the 
Helmholtz  function  for  fluid  water. 


where  the  independent  variables  are  density  (p)  and  temperature  (T). 

Eq.  (1)  is  what  we  mean  when  we  refer  to  the  thermodynamic  surface  for 
water. 

3.  The  Helmholtz  function 

Presented  in  this  section  are  the  equations  that  define  the  Helmholtz 
function  Eq  (1).  The  units  used  for  the  independent  variables  are 
p(g/cm3)  and  T(K),  and  R (J  g_1  K"1)  for  the  specific  gas  constant,  so 
that  A(p,T)  is  given  in  joules/g. 

The  base  function: 


where  a,  £,  and  y are  geometric  constants  11,  133/3  and  7/2;  y is  a 
dimenionless  density,  given  by  y = bp/4,  and  PQ  = .101325  MPa  = 1 atmosphere. 


0) 


• 5Z|±1  + tn  £j*I 


(2) 


0 
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The  parameter  b is  a temperature  dependent  hard  sphere  volume,  and  B is 
related  to  the  second  virial  coefficient. 


where  T = 647.073  K and  the  coefficients  b.  and  B.  are  listed  in  table  1. 
o 11 


bi(cm3  g'1) 

Table  1 
i 

B.j  (cm3  g"1) 

.7478629 

0 

1.1278334 

-.3540782 

1 

- .5944001 

0 

2 

-5.010996 

.007159876 

■ 3 

0 

0 

4 

.63684256 

-.00352836 

5 

0 

The  residual  function: 


+ 


* 


exp 


(5) 
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where  a = 1 cm3/g,  and  the  g. 
data.  The  quantities  6^  and  x. 
respectively,  given  by 


■ 


P - p. 


are  coefficients  determined  by  fits  to 
are  reduced  densities  and  temperatures 
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and  Pi,  T.  are  specified  density  and  temperature  values.  The  k(i)  and 
n(i)  are  integers.  Values  for  the  constants  and  parameters  for  the 
residual  function  are  listed  in  table  2. 


Table  2 


1 

m) 

1(1) 

S(i) 

i 

3c  ( 1 > 

1(1) 

fi  ( i ) 

(J  9’1) 

(J  g’1) 

1 

i 

l 

-.530629685290+03 

19 

5 

4 

-.138025771779+07 

2 

i 

2 

.227449014244+04 

20 

5 

6 

-.251099143690+06 

3 

i 

4 

.787793330207+03 

21 

6 

1 

.465616261156+07 

4 

i 

6 

-.698305273750+02 

22 

6 

2 

-.727527732754+07 

G. 

w 

2 

1 

.178638328754+05 

23 

6 

4 

.417742461483+06 

6 

2 

2 

-.395147315633+05 

24 

6 

6 

.140163582446+07 

7 

2 

4 

.338038842808+05 

25 

7 

1 

-.315552313921+07 

8 

2 

6 

- . 138550502027+05 

26 

7 

2 

.479296663846+07 

9 

3 

1 

-.256374366133+06 

27 

7 

4 

.409126647812+06 

10 

T 

w 

2 

.482125759814+06 

28 

7 

6 

-.136263693884+07 

11 

3 

4 

-.341830169697+06 

39 

9 

1 

.696252208627+06 

12 

•7 

w 

6 

.122231564174+06 

30 

9 

2 

-.108349000964+07 

13 

4 

1 

.117974336556+07 

31 

9 

4 

-.227228274017+06 

14 

4 

2 

-.217348101104+07 

32 

9 

6 

.383654860007+06 

15 

4 

4 

.108299521686+07 

33 

3 

0 

.688332579443+04 

16 

4 

6 

-.254419980640+06 

34 

3 

3 ~ 

.217572455226+05 

17 

5 

1 

-.313777749478+07 

35 

1 

3 

-.266279448298+04 

18 

5 

2 

.529119107577+07 

36 

5 

3 

-.707304180821+05 

i 

MO 

MO 

P-j(g/cm-3) 

T-j  ( K) 

ai 

6i 

9-j  (J  g-1) 

37 

38 

2 

0 

9 

.319 

640. 

34 

2000Q 

-.225 

39 

C 

C 

r\ 

. 319 

640. 

40 

20000 

-1.68 

40 

L 

A 

U 

. 319 

*1  r r 

641.6 

30 

40000 

.055 

*T 

U 

1 . 55 

270. 

1050 

25 

-93.0 
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The  ideal  gas  function: 


where  TR  = T/100  K. 

The  coefficients  C.  are  given  in  table  3. 


Table  3 


1 .197372710180+02 

2 .209662661977+02 

3 -.483429455355+06 

4 .605743189245+05 

5 .225602388500+02 

6 -.987532442000+01 

7 -.431355365532-01 

8 .458155781927-04 

9 -.477549017624-07 

10  .412384608402-10 

11  -.279290527404-13 

12  .144816952031-16 

13  -.564736589529-20 

14  .162004460052-23 

15  -.330362277656-27 

16  .451916066716-31 

17  -.370734122641-35 

18  .137546067535-39 


The  Eqs  1 thru  6 and  the  values  for  the  coefficients  and  parameters 
included  in  this  section  contain  the  complete  thermodynamic  description 
for  water.  Since  A(p,T)  is  everywhere  analytic,  it  is  straightforward 
to  evaluate  appropriate  derivatives  and  to  obtain  closed-form  expressions 


for  all  thermodynamic  properties  over  the  entire  fluid  range.  In  the 
next  section  is  a discussion  of  how  this  is  done. 

4.  Method  of  calculation  and  the  computer  program 

The  thermodynamic  quantities  are  calculated  from  A(p,T)  using  the 
following  relations: 


p 

s 

pRTZ  = p2  3A/3p 

(Pressure) 

(7) 

3P 

3p 

= 

- p + p2 

P 3P2 

(8) 

3P 

3T 

= 

2 32A 
p 3p  3l* 

(9) 

S 

= - 

3A 

3T 

(entropy) 

(10) 

U 

= 

A + TS 

(internal  energy) 

(11) 

H 

= 

U + P/p 

(enthalpy) 

(12) 

Cv 

= - 

T 32A 
1 3T2 

(isochcric  heat  capacity) 

(13) 

G 

= 

A + P/p 

(Gibbs  functions) 

(18) 

cp 

= 

r . T (3P/3TJ2 
Lv  p2  3P/3p 

(isobaric  heat  capacity) 

(19) 

O) 

- 1 

(vcv  fY/2 

(speed  of  sound) 

(20) 

Bn 

l h2P\ 

2RT  V 3p 2]  p=0 

(2nd  virial  coef.) 

(21) 

dBn 

dl 

(22) 

6t 

- | 

I3h\  = 1 _ T 3P/3T 

^P yj  P p'2  3 P/  3p 

(isothermal  J-T  coef.) 

(23) 

5h 

■ ( 

(3T\  = 

^P/H  Cp 

(Joule-Thomson  coef.). 

(24) 
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where  except  as  indicated  otherwise  all  dependent  variables  are  functions 
of  p and  T.  Since  Eqs  7 thru  18  are  linear  in  A(p,T)  each  can  be  represented 
as  the  sum  of  contributions  from  its  three  parts. 

The  computer  print- out  in  appendix  A is  the  FORTRAN  77  program  for 
calculating  thermodynamic  properties  based  on  the  equations  contained  in 
this  report.  All  calculations  related  to  the  same  part  of  A(p,T)  are 
included  in  the  same  sub- routine.  Thus  the  sub- routine  "BASE"  includes 
Eqs  (2-4)  for  the  contribution  from  the  base  function  to  the  Helmholtz 
function,  plus  each  of  the  contributions  therefrom  to  Eqs  7 thru  18. 
Similarly  the  sub- routine  "QQ"  includes  all  such  contributions  from  the 
residual  function,  and  sub-routine  "IDEAL"  contains  the  contributions 
from  the  ideal  gas.  Because  the  detailed  equations  are  easily  "read" 
from  the  respective  sub-routines  they  are  not  reproduced  here. 

The  contributions  to  the  various  thermodynamics  properties  Eqs  7-24 
are  summed  in  sub-routine  "THERM".  Also  included  are  routines  for 
properties  associated  with  liquid-vapor  coexistence.  All  calculations 
are  referred  to  the  liquid  at  the  triple  point,  for  which  state  the 
internal  energy  and  entropy  are  zero.  Though  the  units  for  the  calculations 
are  p(g/cm3),  T(K)  and  joules,  routine  "UNIT"  is  included  which  allows 
use  of  other  units.  Lastly,  the  iteration  sub-routine  "DFIND"  allows 
access  with  independent  variables  P and  T. 

All  equations,  including  parameters  and  constants  listed  in  sections 
2 thru  4 are  included  in  the  program  listing  in  Appendix  A.  Should  any 
inconsistencies  be  found  between  the  text  and  Appendix  A,  it  is  likely 
that  the  errors  are  in  the  text.  Please  let  us  know  if  you  discover 


any. 
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5.  Discussion 

The  data  selected  for  the  derivation  of  the  thermodynamic  surface 
were  primarily  PpT  data.  These  cover  the  range  from  0°  <_  t <_  1000°C  and 
0 < P < 1000  MPa.  Other  kinds  of  data  used  include  values  for  the 
enthalpy  of  the  saturated  liquid  and  for  the  isothermal  compressibility 
of  the  liquid  below  100°C.  No  mathemtaical  constraints  were  used  to 
impose  exact  accord  with  pre-assigned  values.  The  surface  yields  values 
consistent  with  all  input  data  within  reasonable  estimates  of  their 
tolerance.  It  has  been  further  validated  by  extensive  comparison  with 
high  quality  thermodynamic  data,  including  other  than  PpT  that  had  not 
been  used  in  the  derivation  of  the  surface. 

At  very  high  pressure  and  at  very  high  temperatures  thermodynamic 
values  calculated  with  the  surface  are  in  accord  with  those  given  by  the 
(theoretical)  base  function,  so  that  the  surface  should  yield  useful 
estimates  well  beyond  the  range  of  the  data  used  in  its  derivation. 

Thus,  except  very  close  to  the  melting  curve,  the  surface  should  remain 
quantitative  at  pressures  equal  to  2000  MPa  and  at  least  semi-quantitative 
to  4000  MPa,  for  temperatures  in  the  range  250  ± t <_  2500  K for  undissociated 
fluid  water.  Lastly,  because  the  surface  is  analytic,  it  yields  some 
results  at  the  critical  point  and  its  immediate  neighborhood  which  are 
at  variance  with  what  modern  theory  predicts.  Though  we  have  not  found 
significant  departure  from  any  reliable  critical  region  data,  we  caution 
that  the  surface  may  not  conform  to  all  theoretical  expectations  in  the 
region  defined  by 

T = T + 0.5  K 
c — 

p = Pc  (1  ± -2)  , 

where  Tc  and  pc  are  values  at  the  critical  point. 


' 


. 

. 


- 


REFERENCES 


1.  L.  Haar,  J.G.  Gallagher,  and  G.S.  Kell,  in  preparation. 

2.  IUPAC  Commission  on  Atomic  Weights  and  Isotopic  Aburdances,  Pure 
and  Applied  Chem.  52,  2349  (1980). 

3.  E.R.  Cohen  and  B.N.  Taylor,  J.  Phys.,  Chem.  Ref.  Data  2 , 663. 
(1973). 

4.  L.  Haar  and  S.H.  Shenker,  Proc.  Fifth  Symp.  on  Thermophys.  Prop. 
A.S.M.E.  (1970) 

5.  L.  Haar  and  S.H.  Shenker,  J.  Chem.  Phys.  ^5  4951  (1971). 

6.  J.R.  Wardell,  G.  Wilmot,  L.  Haar  and  M.  Klein,  Proc.  12th  JANNAF 
Combustion  Meeting,  Chem.  Prop.  Info.  Agency,  publ . 273  (1975). 

7.  L.  Haar,  G.  Powell,  M.  Klein,  J.  Wardell,  and  G.  Wilmot,  Proc.  13th 
JANNAF  Combustion  Meeting,  Chem.  Prop.  Info.  Agency,  publ.  281 
(1976). 

8.  G.  Powell,  G.  Wilmot,  L.  Haar  and  M.  Klein,  Interior  Ballistics  of 
Guns,  M.  Summerfield,  Ed.,  Part  IV,  P.  307,  Am  Inst.  Aero,  and 
Astro.,  New  York  (1979). 

9.  M.  Klein,  L.  Haar,  G.  Powell,  and  G.  Wilmot,  15th  JANNAF  Combustion 
Meeting,  Chem.  Prop.  Info.  Agency,  Publ.  297  (1978). 

10.  F.  Kohler  and  L.  Haar,  accepted  for  publ.  in  J.  Chem.  Phys. 

11.  L.  Haar,  J.  Gallagher,  in  preparation 

12.  L.  Haar,  J.  Gallagher,  and  G.  Kell,  Proc.,  9th  Int.  Conf.  on  Prop. 
Steam,  Munich  (1979),  J.  Straub  and  K.  Scheffler,  Eds.,  P.  69, 
Pergamon  Press,  Oxford  (1980). 

13.  H.W.  Woolley,  Proc.,  9th  Int.  Conf.  on  Prop.  Steam,  Munich  (1979), 
J.  Straub  and  K.  Scheffler,  Eds.,  P.  166,  Pergamon  Press,  Oxford 
(1980). 


' 


■ 


' 


It 


APPENDIX  A 


This  appendix  contains  the  FORTRAN  77  programs  which  yield 
values  for  thermodynamic  properties  for  liquid  and  gaseous  states 
for  water.  The  independent  variables  are  temperature  (T)  and 
density  (D).  The  programs  should  yield  useful  results  in  the  range 

250  <_  T(K)  < 2500 

10“8  <_  D(g/cm3)  <_  D(Pmax), 

where  Pmax  is  the  lesser  of  the  pressure  of  melting  ice  or  4000  MPa. 
For  values  of  density  less  than  or  equal,  to  10" 8 g/cm3,  the  program 
sets  the  density  of  10“ 8 g/cm3. 
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Outlined  below  are  the  procedures  for  use  of  the  subroutines  which 
follow.  In.this  outline,  the  symbol  used  tfor  pressure  is  "P’  , for 
density  is  "B” , and  for  temperature  is  ”T".  Calculations  are  made  in  the 
internal  units  of  the  program  which  are  MPa,  g/cm3  and  deg  K respectively. 

See  III  below  for  conversions  to  other  units.  Examples  of  the  use  of  these 
routines  are  given  in  three  main  programs  in  Appendix  B,  with  samples  of 
the  output  which  can  be  used  for  checking  the  implementation  of  the  programs. 

I.  Jor  Density  and  Temperature  as  independent  variables,  the  sequence 

CALL  BB(T) 
call  base (d , t ) 

CALL  Q0(T ,D) 

will  place  quantities  in  COMMON  from  which 
P s D^GASCON^T^ZB  + Q0 
dP/dD  = GASCON*.?*  (ZB+T*DZB)  + Q5 
To  obtain  the  other  thermodynamic  properties: 

CALL  THERM (D,T) 

will  olace  in  COMMON  "AB“  (=A/RT),  "GB"  (=G/RT ) r "SB”  (=S/R), 

"UB"  (=U/RT),  "HB“  (=H/RT),  ”C7B"  (=Cv/R),  ‘*CPD'  £=Cp/R), 

DPET'  ( =dP/dT , MPa/K),  "DVD?"  (=d7/dT,  cm3/gl) , ‘ CJTT" 

(isothermal  joule-thompson  coeff,  cm3/g)  and  ’CJTH"  (adiabatic 
joule-thompson  coeff,  K/MPa).  These  properties  in  the  desired 
units  can  then  be  obtained  by  R or  RT  as  appropriate  in  the 
desired  units  . 

II.  Pressure  and  Temperature  as  independent  variables 

CALL  DJIND( DOUT , F , LGUESS  ,T,DPDD) 

This  will  return  the  density  corresponding  to  the  input  P and  T as 
DOUT  from  which  the  procedures  in  I can  be  used  for  further  calcula- 
tions. This  routine  requires  an  initial  approximate  guess  for  the 
density  for  use  to  begin  the  Newton  iteration.  dP/dD  is  also  returned 
as  a byproduct. 

III.  Use  of  other  systems  of  units. 

A group  of  three  subroutines  is  included  which  facilitate  the  use 
of  a system  of  units  chosen  by  the  user: 

CALL  UNIT 

will  query  the  user  as  to  which  units  are  desired  for  temperature, 
pressure,  density  and  energy,  and  set  up  the  necessary  parameters  for 
converting  from  the  desired  units  to  the  internal  units  and  baci 
again  for  the  output.  The  names  of  the  units  in  alphanumeric 
characters  are  also  placed  in  COMMON  for  use  in  laoeling. 

FUNCTION  TTT ( Tex ternal ) , and 
JUNCTION  TTI ( Tin te rnal ) 

will  convert  from  external  T to  deg  £,  and  from  deg  K to  the  external 
T respectively. 


O O O o O 


A3 


BLOCK  DATA 

implicit  double  precision (a-h, o-z ) 


real  p,q 
COMMON  /ACONST/ 

COMMON  /NCONST/ 

COMMON  /ELLCON/ 

COMMON  /BCONST/ 

COMMON  /ADDCON / 

THIS  BLOCKDATA  SUBROUTINE 
USED  IN  THE  REST  OF  THE 


VM, GASCON, TZ,AA, ZB, DZB ,YB ,UREF ,SREF 
G(40)  ,11(40 ),JJ( 40) ,NC 
G1,G2,GF,B1 ,B2,B1T ,B2T ,B1TT ,B2TT 
P ( 10  ) , 0 ( 10 ) 

ATZ (4 ) , ADZ (4 ) , AAT ( 4 ) ,AAD(4) 

SUPPLIES  MOST  OF  THE  FIXED  PARAMETERS 
ROUTINES.  P IS  THE  b (i ) of  TABLE  I, 


I , G1,G2  AND  GF 
G , 1 1 , JJ  ARE  THE 


ARE  THE  ALPHA, 
g ( i ) , k(i)  AND 


BETA 


Q IS  THE  B ( i ) OF  TABLE 
AND  GAMMA  OF  EC  2,  AND 
1 ( i ) of  EQ  5. 

DATA  ATZ/2*64.D1,641.6D0,27.D1/,ADZ/3*.319D0,1.55D0/,AAT/2*2.D4 
1 ,4.D4,25 . D0/ , AAD/34 . D0 ,4 .Dl ,3 . Dl , 1 .05D3/ 

DATA  WM/18.0152D0/, GASCON/. 461522D0/ , TZ/647 .073D0/ , AA/1 .D0/ , NC /36/ 
DATA  UREF ,S REF/ -43 28 .455039D0 ,7 .6160802D0/ 

DATA  Gl,G2,GF/ll.d0 ,44 .333333333333d0 , 3 . 5d0/ 

DATA  P/.747862S ,-.3540782,2*0. , .007159876  ,0 . ,-.003528426,3*0./ 

DATA  Q/l . 1278334,0 . ,-.5944001,-5 .010996,0. , .63684256  ,4*0  ./ 

DATA  &/- . 52062S68529023D3 , .22744901424408D4 , .78779333020687D3 

1, -. 69830 5273749S4D2, . 17863832875422D5 ,-.39514731563338D5 

2,  .33803884280753D5  ,- . 13855050202703D5 , - .25637436613260D6 

3. . 4821257 598 1415D6 ,- . 34183016969660D6  , .122231564 17448 D6 

4. . 1 1797433655832D7,-. 217348 101 10373D7,  . 1082S952168620D7 

5 , -. 25441998 064 049D6 ,-. 3 1377774 947767 D7  , .52911910757704D7 

6 , - . 13802577 177877D7 , - . 25 10991436900 1D6 , .465618261 15608D7 
7, 72752773275387D7, .41774246148 294D6, . 1401 6358244614D7 

6 ,- .31 555231392127D7 , . 47929666384584D7 , . 40912664 781209D6 
9 ,- . 13626 36S388386D7  , .69625220862664D6 ,- .10 834900096447D7 

A , - . 2272282740 1688D6 , .38365486000660D6 , .63833257944332D4 

B ,  .21757245522644D5 ,-. 2662794482977 0D4,-. 707 30 4180820 74D 5 

C , -. 225D0 , -1 .68D0 , .055D0 ,-93.0D0/ 

DATA  11/4*0 ,4*1 ,4*2, 4*3 ,4*4, 4* 5, 4*6, 4*8, 2*2 ,0,4, 3*2, 4/ 

DATA  JJ/ 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7, 2, 3, 5, 7 
1,2, 3, 5, 7, 1,3*4, 0,2, 0,0/ 

END 
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SUBROUTINE  BB (T ) 

THIS  SUBROUTINE  CALCULATES  THE  B'S  OF  EQS  . 3,4  USING  COEFFICIENTS 
FROM  BLOCXDATA,  CALCULATING  ALSO  THE  FIRST  AND  SECOND  DERIVS 
W.R.TO  TEMP.  THE  B'S  CALCULATED  HERE  ARE  IN  CM3/G . 
implicit  double  precision (a-h, o-z) 
real 

COMMON  /ELLCON/  G1  ,G2 ,GF ,B1 ,E2  ,E1T ,E2T ,B1TT ,B2TT 
COMMON  /ACONST/  WM , GASCON , TZ , AA , Z ,DZ ,T ,UREF ,SREF 
COMMON  /BCONST/  P(10),Q(10) 

DIMENSION  V(10) 

V(l)=l. 

DO  2 1=2,10 
2 V (I )=VII-1)*TZ/T 

E 1=P( 1 )+F ( 2)*AL0G ( 1 . /V ( 2 ) ) 

52=Q(1) 

E1T=P (2)*V(2) /TZ 
E2T=0 . 

B1TT=0  . 

B2TT=0 . 

DO  4 1=3,10 
B1=B1+P ( I )*V(I-1) 

B2=B2+g(I )*V(I-1) 

31T=B IT- ( 1-2 )*P ( I )*V ( 1-1 ) /T 
B2T=B2T-( I-2)*0(I )*V(I-1)/T 
B1TT=B1TT +P ( I ) * ( I -2 )**2*V ( 1-1 ) /T/T 
4 E2TT=B2TT+Q(I )*(I-2 )**2*V ( 1-1 ) /T/T 
B1TT=E1TT-E1T/T 
B2TT=B2TT-B2T/T 
RETURN 
END 
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SUBROUTINE  BASE (D  ,T ) 

C This  function  calculates  Z [=Pbase/( DRT)  ] , 

C and  also  Abase  ,Gbase , Sbase ,Ubase ,H base ,CVbase  , AND  1/(DRT)  * DP/DT 
C for  the  base  fct  (called  DPDTB) 

C The  AB,GB ,SB ,U3,HB  and  CVB  are  calculated  in  dimensionless  units: 

C AB/RT,  GB/RT,  SB/R,  etc. 

IMPLICIT  REALMS  (A-H.O-Z) 

COMMON  /ELLCON/  G1 ,G2 f GF ,Bl  ,B2  ,E IT  ,B2T  ,B1TT  ,B2TT 
G1,G2  AND  GF  ARE  THE  ALPHA,  BETA  AND  GAMMA  OF  EQ  2,  WHICH  ARE 
SUPPLIED  BY  THE  BLOCKDATA  ROUTINE.  B1  AND  B2  ARE  THE  ’‘EXCLUDED 
VOLUME'  AND  “2ND  VIRIAL”  (EQS  3 AND  4)  SUPPLIED  EY  THE  SUBROUTINE 
BB (T ) , WHICH  ALSO  SUPPLIES  THE  1ST  AND  2ND  DERIVATIVES  WITH 
RESPECT  TO  T ( E1T,B2T ,B1TT ,E2TT  ) . 

COMMON  /BASEF/  AB ,GB , SB , UB , HB  , CVB  , DPDTB 
COMMON  /ACONST/  WM  , GASCON , TZ, A ,Z ,DZ ,Y ,UREF ,SREF 
Y= . 25^31*0 
X=1 .-Y 

Z0=(1 .+G1*Y+G2*Y*Y)/X**3 
Z=Z0+4.*Y*(32/B1-GF) 

DZ0=(G1+2.*G2*Y)/X**3  + 3.*(  l.+G  1*Y+G2*Y*Y)  /X**4 
DZ=DZ0+4.*(E2/E1-GF) 

AB  = -DLOG ( X ) - ( G2-1 . )/X+28 . 16666667D0/X/X+4 .^Y* ( B2/B1-GF ) 

1 +15 . 166666667D0  + DLOG ( D*T*4 .55483D0 ) 

G3  = AB  + Z 
BB2TT=T*T*E2TT 

UB=  -T*B1T*(Z-1.-D*B2)/B1-D*T*B2T 
HE=Z+UB 

C VB=2 .*UB+ ( Z0-1  .)*( (T*B1T/B1)**2-T*T*B1TT/B1) 

1 - D*(BB2TT  - GF^BITT^T^T ) -( T*B1T/B1 )**2*Y*DZ0 
DPDTB  = Z/T  + D* ( D Z* B IT /4.+ 323-32/31*3  IT ) 

SB  = UB  - AB 
RETURN 
END 
C 
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SUBROUTI NS  QQ(T  ,D ) 

THIS  ROUTINE  CALCULATES , FOR  A GIVEN  T(K ) AND  D(G/CM3),  THE  RESIDUAL 
CONTRIBUTIONS  TO:  PRESSURE  (Q),  HELMHOLTZ  FCT  (AR)  , DP/DRHO  (Q5), 
AND  ALSO  TO  THE  GIBBS  FUNCTION,  ENTROPY , INTERNAL  ENERGY,  ENTHALPY, 
ISOCEORIC  HEAT  CAFACITY  AND  DPDT . (EQ  5) 

TERMS  37  THRU  39  ARE  THE  ADDITIONAL  TERMS  AFFECTING  ONLY  THE 
IMMEDIATE  VICINITY  OF  THE  CRITICAL  POINT,  AND  TERM  40  IS  THE 
ADDITIONAL  TERM  IMPROVING  THE  LOV  T,  HIGH  P REGION. 

IMPLICIT  REAL*8(A-H,0-Z) 

COMMON  /RESF/  AR,GR ,SR ,UR ,HR, CVR ,DPDTR 
COMMON  /QQQQ/  Q,Q5 

DIMENSION  QR(  11)  ,QT(10)  ,0.ZR(9)  ,QZT(9) 

EQUIVALENCE  (QR(3) ,QZR(l) ) , (QT(2) ,QZT(1) ) 

COMMON  /N CONST / G (40  ) , 1 1 ( 40 ) , J J( 40 ) , N 

COMMON  /ACONST/  WM , GASCON , TZ , AA , Z , DZ ,Y ,UREF , SREF 

COMMON  /ADDCON / ATZ (4 ) ,ADZ (4 ) , AAT (4 ) , AAD (4 ) 

RT=GASCON*T 
QR( 1 ) =0 . 

QS=0. 

Q=0 ,D0 
AR=0.D0 
DADT=0  . 

C VR=0  . 

DPDTR=0 . 

E=DEXP (-AA^D ) 

O10=Da!cD*  E 
Q20=l .D0-E 
QR(2)=Q10 
V=TZ/T 
QT ( 1 ) =T/TZ 
DO  4 1-2,10 
QR(I+1)=QR(I )*Q20 
4 QT(I)=QT(I-1)*V 

DO  10  1=1, N 
K=I I ( I )+l 
L=JJ ( I ) 

Z z = K 

QP=G ( I )* AA*QZ R ( K-l  )*OZT (L  ) 

Q=Q+QP 

C5  = 05  + AA*(2./D-AA*(1.-E*(K-1)/Q20)  )*QP 
AR=AR+G( I )*QZR(K)*QZT(L) /Q10/ZZ/RT 
DFDT=C20**K* ( 1-L ) *QZT ( L + 1 )/TZ /K 
D2F=L*DFDT 

DPT=DFDT*Q10*AA*K/Q20 
DADT=DADT+G ( I )*DFDT 
DPDTR=DPDTR+G ( I )*DPT 
10  CVR=CVR-G(I )*D2F/GASC0N 
CP=0. 

Q2A=0  . 
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DO  20  J=37 ,40 
If(G(J).EQ.0.D0)  GO  TO  20 
K=I  I ( J ) 

KM= JJ { J ) 

DDZ  = AD  Z (J-36 ) 

DEL  - D/DDZ  ~ 1 

IE (DABS (DEL) .LT.l.D-10)  DEL=1.D-10 
DD  = DEL-DEL 
EX1  = -AAD( J-36)*DEL**K 
DEX=DEXP(EX1 )*DEL**KM 
ATT  = AAT ( J-36 } 

TX  = ATZ( J-36  ' 

TAU  = T/TX-1 . 

EX2  = -ATT^TAU-TAU 
TEX  = DEXP ( EX2 ) 

CIO  = TVS’ TTY 

QM  = KM/ DEL  - K*AAD (J-36 j^DEL** ( K-l ) 
FCT=QM*D**2*Q10/DDZ 

Q5T  = FCT*(2./D+QM/DDZ) -(D/DDZ  )**2*Q10* ( KM/DEL /DEL- 
1 K* (K-l )-AAD ( J-36  j^DEL** ( K-2 ) ) 

05  = Q5  + QST^G ( J ) 

CP  = OP  + G ( J )*ECT 

DADT  = DADT  - 2 .*G ( J )*ATT*TAU*O10/TX 
DPDTR  = DPDTR  - 2 .*G ( J)*ATT*TAU*FCT/TX 
02A  = Q2A  + T*G(J)*(4.*ATT*EX2+2.*ATT)*O10/TX/TX 
AR  = AR  + 010-G ( J ) /RT 
20  CONTINUE 

SR=-DADT/GASCON 

UR=AR+SR 

C VR=C VR+  Q2 A/GAS  CON 

Q=0+QP 

RETURN 

END 


C 
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SUBROUTINE  DFIND (DOUT,P , D,T,DPD) 

ROUTINE  TO  FIND  DENSITY  CORRESPONDING  TO  INPUT  PRESSURE  P(MPA),  AND 
TEMPERATURE  T(K),  USING  INITIAL  GUESS  DENSITY  D(G/CM3).  THE  OUTPUT 
DENSITY  IS  IN  G/CM3 , ALSO,  AS  A BYPRODUCT,  DP/DRHO  IS  CALCULATED 
( "DPD" , MPA  CM3/G) 

IMPLICIT  REAL*8 ( A-H ,0-Z ) 

COMMON  /QCQQ/  Q0,Q 5 

COMMON  /ACONST/  WM, GASCON, TZ,AA,Z,DZ,Y, UR EF ,SREF 
DD=D 

RT=GASCON5!‘T 
IF( DD . LE .0 . ) DD=1 .D-8 
IF ( DD . GT . 1 . 9 ) DD=1 . 9 
L=0 

9 L=L+1 

11  IF(DD.LE.0.)  DD=l.D-8 
IF( DD .GT . 1 .9 ) DD=1 .9 
CALL  BASE  (DD , T ) 

CALL  QQ(T  ,DD) 

PP  = RT*DD*Z  + Q0 
DPD=RT*(Z+Y*DZ)+Q5 
C THE  FOLLOWING  3 LINES  CHECK  FOR  NEGATIVE  DP/DRHO,  AND  IF  SO  ASSUME 
C GUESS  TO  BE  IN  2-PHASE  REGION,  AND  CORRECT  GUESS  ACCORDINGLY. 

IF(  DPD  .GT  .0  ,D0 ) GO  TO  13 
IF(D.GE. .2967D0 ) DD=DD*1.02D0 
IF(D.LT.  .2967D0)  DD=DD*.98D0 

IF (L.LE. 10 ) GO  TO  9 * 

13  DPDX=DPDsicl . 1 

IF ( DPDX . LT . . 1 ) DPDX= . 1 
DP=DABS ( 1 .-PP/P ) 

IF(DP .LT.l.D-e)  GO  TO  20 
IF ( D . GT . .3  .AND.  DP.LT.l.D-7)  GO  TO  20 
I F ( D . GT . .7  .AND.  DP.LT.l.D-6)  GO  TO  20 
X= ( P-PP ) /DPDX 

IF(DABS(X)  .GT  . .1)  X=X* . 1/DABS (X ) 

DD=DD+X 

IF ( DD . LE . 0 . ) DD=l.D-8 
19  IF ( L . LE. 30 ) GO  TO  9 
20  CONTINUE 
DOUT=DD 
RETURN 
END 
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SUBROUTINE  THERM(D,T) 

IMPLICIT  REAL*8 (A-H  0”Z ) 

COMMON  /ACONST/  WM , GASCON , TZ , AA , ZB ,DZB ,Y ,UREF ,SREF 
COMMON  /QOQQ/  QP , QDP 

COMMON  /EASEF/AB ,GB , SB ,UE ,HB , C VB  ,DPDTB 
COMMON  /RESF/AR ,GR  , SR, UR , HR , CVR , DPDTR 
COMMON  /IDF/  AI  ,GI ,SI,UI,HI,CVI,CPI 

COMMON  /FCTS/AD ,GD , SD,UD ,HD  , CVD , CPD ,DPDT  ,DVDT , C JTT , CJTH 
THIS  SUBROUTINE  CALCULATES  THE  THERMODYNAMIC  FUNCTIONS  IN 
DIMENSIONLESS  UNITS  (AD=A/RT , GD=G/RT,  SD=S/R,  UD=U/RT, 

HD=H/RT,  CVD=C V/R , AND  CPD=CP/R) 

CALL  IDEAL (T) 

RT  = GASCON*T 
Z = ZB  + QP/RT/D 
DPDD  = RT* ( ZB+Y*DZB ) + QDP 
AD  = AE  + AR  + AI  - UREF/T  + 5REF 
GD  = AD  + Z 

UB  = UB  + UR  + UI  - UREF/T 
DPDT  = RT*D*DFDTB  + DPDTR 
CVD  = CVB  + CVR  + CVI 

CPD  = CVD  + T*DPDT**2/(D*D*DPDD*GASC0N) 

HD  = UD  + Z 

SD  = SB  + SR  + SI  - SREF 
DVDT=DPDT /DPDD /D/D 
CJTT=1./D-T*DVDT 
C JTH^-CJTT/CPD/GAS  CON 
RETURN 
END 
C 

FUNCTION  PS ( T ) 

C This  function  calculates  an  approximation  to  the  vapor  pressure,  PS, 

C as  a function  of  the  input  temperature.  The  vapor  pressure 
C calculated  agrees  with  the  vapor  pressure  predicted  by  the*  surface 

C to  within  .02?  to  within  a degree  or  so  of  the  critical  temperature, 

C and  can  serve  as  an  initial  guess  for  further  refinement  by 
C imposing  the  condition  that  Gl=Gv.  t 

IMPLICIT  REAL*6  (A-H.O-Z) 

DIMENSION  A (8 ) /-7 .8889166D0 ,2 . 5514255D0 ,-6 .716169D0 
1 ,33 .23949 5D0, -105 .38479D0 , 174 .35319D0 ,-148 .39348D0 
2 , 48 .631602D0/ 

IF ( T . GT . 214  .D2 } GQ  TO  2 

PL=6 . 3573118D0-8858 .843D0/T+607 . 56335D0*T** (-.6) 

PS=.1*D£XP(PL) 

RETURN 

2 V=T/ 647 . 25D0 
W=DABS ( 1 .D0-V ! 

B=0 .D0 
DO  4 1=1,8 
Z = I 

4 B=B+A(I)*W**((Z+l.)/2.) 

Q=B/V 

PS=22 . 093D0*BEXP ( Q ) 

RETURN 

END 


A 


C 

FUNCTION  TSAT(P) 

C This  function  calculates  the  saturation  temperature  for  a given 
C pressure,  by  an  iterative  process  using  PSAT  and  TDPSDT. 

REAL*8  P , PS  , TO,  TSAT, TDPSDT 
TSAT=0 . 

IF (P .GT .22.05)  RETURN 
K=0 

PL=2 . 3025£5*DLOG ( P ) 

C PL=LOGE(10)*LOGE(P)  TO  CONVERT  EQUATION  FROM  BARS  TO  MPa 

TG=372 .8  3+?!* ( 27 .7589+PL* (2  .3819+PL* ( .24834 +PL* .0193855 ) ) ) 

1 IF (TG  .LT  . 273 . 15 ) TG=273.15 
IF( TG .GT .647 . ) TG=647. 

I F ( K . LT . 8 ) GO  TO  2 
WRITE (6, 3 ) K , P ,PP , TG 

3 FORMAT () 

GO  TO  8 

2 K=K  + 1 
PP=PS ( TG ) 

DP= TDPSDT (TG ) 

IF(ABS(1.-PP/P) .LT. .00001)  GO  TO  8 
TG  = TG* ( 1 . + ( P-PP ) /DP ) 

GO  TO  1 
8 TSAT=TG 
RETURN 
END 
C 

FUNCTION  TDPSDT (T ) 

C This  function  calculates  T*(dPs/dT),  and  is  used  by  the  function  TSAT 
IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  A( 8 ) /-7 .888S166D0 , 2 . 5514255D0 , -6. 716169D0 
1 ,33 .239495D0 ,-105 .38479D0 ,174.35319D0 ,-148 . 39348D0 
2,48 .631602D0/ 

V =T/ 647. 25 
W=1  .-V 
B=0  . 

C=0  . 

DO  4 1=1,8 
Z=I 

Y=A(I }*¥**( (Z+l.J/2. ) 

C=C+Y/W*( .5-.5*Z-l./V) 

4 B=B+Y 
Q=B/V 

TDPSDT =2 2 . 093*DEXP (Q )*C 

RETURN 

END 
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SUBROUTINE  IDEAL(T) 

C THIS  SUBROUTINE  CALCULATES  THE  THERMODYNAMIC  PROPERTIES  FOR 
C WATER  IN  THE  IDEAL  GAS  STATE  OF  H .W  . WOOLLEY 
IMPLICIT  REALMS  (A-H.O-Z) 

COMMON  /IDF/  Al  ,GI ,SI,UI,HI,CVI,CPI 

DIMENSION  C( 18)/.19720271018D2f . 209662681977D2 483429455355D0 

1 . . 605743 189245D1 ,22 .560 2388 5D0  ,-9 .87532442D0 ,- .43135536513D1 

2. . 458155781D2 , - .47754901883D-1 , . 41238460633D-2 27929052352D-3 

3 . .  1448 159526 1D -4 , - . 56473658748D-6 , .16200446D-7,-.3303822796D-9 

4. . 451516067368D-11 . 370734122706D-13 , . 137546068238D-15/ 

C NOTE  THAT  THE  TEMPERATURES  ARE  SCALED  BY  A FACTOR  OF  100  HERE  SO  THAT 
C THE  EXPONENT  OF  THE  COEFFICIENT  OF  THE  LAST  TERM  WILL  BE  WITHIN 
C RANGE  FOR  MOST  COMPUTERS. 

TT=T/1  .D2 
TL=DLOG( TT ) 

GI=-(C(1)/TT+C(2) )*TL 
HI=(C(2)+C(l)*(l.D0-TL)/TT) 

CPI =C ( 2 ) -C ( 1 ) /TT 
DO  8 1=3,16 
GI=GI-C(I )*TT**(I-6'> 

HI  = HI+C(I  )*(I-6)*TT**(I-6) 

8 CPI=CPI+C(I)*(I-6)*(I-5)*TT**(I-6) 

AI=GI-1. 

UI=HI -1 . 

CVI=CPI-1. 

SI=UI-AI 

RETURN 

END 

C 

SUBROUTINE  SECVIR ( T , VIR  ' 

C THIS  SUBROUTINE  CALCULATES  THE  SECOND  VIRIAL  IN  CM3/G 
C AT  TEMPERATURE  T IN  K. 

IMPLICIT  REAL*8 (A-H.O-Z) 

COMMON  /N  CON  ST / G (40 } , 1 1 (40  ), J J ( 40 ) , NC 

COMMON  /ELLCON/  G 1 ,G2 , GF , BB1 , BB2 ,B IT ,B2T ,B1 TT , B2TT 

COMMON  /OQCQ/  02,05 

COMMON  /ACONST/  WM .GASCON , TC ,AA , Z , DZ ,Y  ,UREF ,SREF 
CALL  BB( T ) 

V = T C / T 
V IR=BB2 
DO  3 J=1  ,NC 

IF ( 1 1 ( J ) . ME  .0  } GO  TO  3 
L=JJ ( J ) 

VIR=VIR+G( J)*V**(L-1  '/T /GASCON 
3 CONTINUE 
RETURN 
END 
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SUBROUTINE  CORR (T ,P ,DL,DV,DELG  ) 

C SUBROUTINE  CORR  WILL  CALCULATE,  FOR  AN  INPUT  T AND  P AT  OR 
C VAPOR  PRESSURE,  THE  CORRESPONDING  LIQUID  AND  VAPOR  DENSIT 
C DELG  = ( GL-GV ) /RT  FOR  USE  IN  CALCULATING  THE  CORRECTION  T 
C PRESSURE  FOR  DELG  = 0. 

IMPLICIT  DOUBLE  PRECIS  ION (A-H , O-Z ) 

COMMON  /QQQQ/  Q00,Q11 

COMMON  /ACONST/  WM , GASCON , TZ , AA , ZB ,D ZB ,TB ,UREF , SREF 
COMMON  /FCTS / AD , GD , SD , UD , HD , C VD , CPD , DPDT , D VDT , C JTT , CJ 
DLIQ=DL 

I F(  DL  . LE  .0  . ) DLIQ  = 1 . 11-.0004*T 
CALL  BB( T ) 

RT=GASCON*T 

CALL  DFI ND (DL  ,P  ,DLIQ ,T , DQ) 

CALL  THERM (DL.T) 

GL=GD 

DVAP=DV 

IF( DV  .LE  .0  . ) DVAP=P/GAS CON /T 
CALL  DFI ND(DV,P,DVAP,T,DQ) 

IF(DV.LT.5.D-7)  DV=5.D-7 
CALL  THERM (DV  ,T) 

GY=GD 

DELG  - GL-GV 
RETURN 
END 
CC 

SUBROUTINE  PCORR ( T ,P  ,DL ,DV ) 

C SUBROUTINE  PCORR  WILL  CALCULATE  THE  VAPOR  PRESSURE  P AND  T 
C VAPOR  DENSITIES  CORRESPONDING  TO  THE  INPUT  T,  CORRECTED  S 
C GL-GV=0 . THE  FUNCTION  PS  IS  REQUIRED  WHICH  WILL  GIVE  A R 
C GOOD  APPROXIMATION  TO  THE  VAPOR  PRESSURE  TO  BE  USED  AS  TH 
C POINT  FOR  THE  ITERATION. 

IMPLICIT  DOUBLE  PRECISION  (A-H, O-Z) 

COMMON  /ACONST/  WM  , GASCON , TZ , AA , ZB ,DZB  ,YB ,URFF , SREF 
P = PS(T) 

2 CALL  CORR (T,P,BL,DV,DELG) 

DP=DELG*GASCON*T/(l./DV-l./DL) 

P = P+DP 

IF ( DABS ( DELG ) .LT.l.D-4)  RETURN 
GOTO  2 
END 
CC 

SUBROUTINE  TCORR( T , P , DL , DV ) 

SUBROUTINE  TCORR  IS  SIMILAR  TO  ''PCORR"  EXCEPT  THAT  THE  TEM 
CORRESPONDING  TO  THE  INPUT  VAPOR  PRESSURE  IS  FOUND.  FUNCT 
ARE  TSAT  AND  TDPSDT  WHICH  GIVE  AN  APPROXIMATION  TO  T ( SAT ) 
T*DP(SAT)/DT . 

IMPLICIT  DOUBLE  PRECISION  (A-H, O-Z) 

COMMON  /ACONST/  WM  , GASCON , TZ , AA , ZB ,DZE , IB ,UREF , SREF 
T = TSAT(P) 

2 CALL  CORR(T,P,DL,DV  ,DELG) 

DP= DELG5*  GAS  CON^T/tl./DV-l./DL) 

T = T* (1 . -DP/TDPSDT (T } ) 

I F ( DA BS ( DELG ) . LT . 1 .D-4 ) RETURN 
GO  TO  2 
END 


c 


c The  following  3 subroutines  form  a package  allowing  the  us 
c operate  in  a s/stem  of  units  of  his  choice. 

SUBROUTINE  UNIT 

C THIS  SUBROUTINE  QUERIES  THE  USER  AS  TO  HIS  CHOICE  OF  UNITS 
C INTERNAL  PARAMETERS  APPROPRIATELY.  THE  INTERNAL  UNITS  OF 
C TEMPERATURES  IN  K , DENSITIES  IN  G/CM3 , ALL  OTHER  QUANTITI 
C IN  DIMENSIONLESS  UNITS  AND  DIMENSIONED  AT  OUTPUT  TIME. 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  NT  , ND , NP , NH, NNT , NND , NNP , NNH 
COMMON  /UNITS/  IT , ID , IP , IH , NT , ND  . NP , NH ,FT  ,FD ,FP , FH 
DIMENSION  FFD ( 4 ) , FFP  ( 5 ) , FFH (6) , N NT ( 4 ) , NND (4  ) , NNP ( 5)rNN 
DATA  FFD/l.D-3, 1.D0, .0180152D0, .016018D0/ 

DATA  FFF/1.D0 , 10.D0 ,9 .869232667D0 , 145 . 038D0 ,10 .1971D0/ 
DATA  FFH/2*1.D0,18.0152D0, . 23884590D0 ,4 .30285666D0 , .42 
DATA  NNT/1HK , IH C , 1HR  , 1HF/ 

DATA  NND/6Hkg/m3  ,6Hg/cm3  ,6Hmol/L  ,6Hlb/ft3/ 

DATA  NNP/6H  MPa  ,6H  Bar  , 6H  Atm  ,6H  PS  I ,6Hkg/cm2/ 
DATA  NNH/6hkJ/kft  ,6H  J/g  f6HJ/'mol  ,6Hcal/g  ,7Hcal/mol 
DATA  Al, A2 , A3 .A4/8HTEMPERAT ,7HDEN5ITY 
1 ,8HPRESSURE , 8HENERGY  / 

WRITE ( 6,11)  Al 

30  WRITE (6,12) 

READ  ( 5 ,10*fEND=99}  IT 
IF( IT . EQ  .0 ) STOP 
IF( IT .GT . 4 ) GOTO  30 
NT=NNT(IT ) 

WRI TE ( 6 , 11 ) A2 

31  WRITE (6,13) 

READ ( 5 ,10 , END=99 ) ID 

IF ( ID . GT . 4 .OR.  ID.LT.l)  GOTO  31 

ND=NND( I D ) 

FD=FFD ( ID ) 

WRITE(6, 11 ) A3 

32  WRITE (6, 14) 

READ ( 5 ,10 ,END=S9)  I? 

IF ( IP  .GT  .5  .OR.  IP.LT.l)  GOTO  32 
N P=  N N F ( I P ) 

FP=FFF (IP) 

WRITE (6,11)  A 4 

33  WRITE(6,15) 

READ( 5,10  , END =99)  IH 

IF ( IH  .GT  . 6 .OR.  IH.LT.l)  GOTO  33 

NH=NNH (IH ) 

FH=FFH ( I H ) 


99 

10 

11 

12 

13 

14 

15 


RETURN 
STOP 
FORMAT ( ) 
FORMAT ( ' 
FORMAT ( ' 
FORMAT ( ' 
FORMAT ( ' 
FORMAT ( ' 
1RIES/MOL 
END 


ENTER  UNITS 
CHOOSE  FROM 
CHOOSE  FROM 
CHOOSE  FROM 
CHOOSE  FROM 
6=BTU/LE ') 


CHOSEN  FOR  ' , A8 ) 

1=DEG  K,  2= DEG  C,  3=DEG  R,  4=DEG 
1 =EG/M3,  2=G/CM3 , 3=M0L/L , 4=LE/F 
1 = MPA , 2=BAR , 3=ATM  , 4=PSIA,  5=KG 
1=KJ/KG,  2=J/G,  3=J /MOL , 4=CALORI 


A'* 


FUNCTION  TTT(T) 

C FUNCTION  TO  CONVERT  INPUT  TEMPERATURES  IN  EXTERNAL  UNITS  T 
DOUBLE  PRECISION  T , TTT ,FT  »FD , FP , FH , NT , ND , NP ,NH 
COMMON  /UNITS/  IT , ID  , IP , IH , NT f ND , NP , NH ,FT ,FD ,FP , FH 
GO  TO  (1,2, 3, 4) ,IT 

1 TTT=T 
FT-1 . 

RETURN 

2 TTT=T -*-273 . 15D0 
FT=1 . 

RETURN 

3 TTT=T/1.8D0 

FT=  *55555 555 555 56D0 
RETURN 

4 TTT=( T+459 .67D0 )/l .8D0 
FT= . 5555555555556D0 
RETURN 

END 

FUNCTION  TTI(T) 

C FUNCTION  TO  CONVERT  INTERNAL  TEMPERATURES  IN  DEG  K TO  EXTE 
DOUBLE  PRECISION  T ,TTI , FT , FD ,FP ,FH , NT  ,ND , NP ,NH 
COMMON  /UNITS/  IT , ID , IP , IH , NT , ND , NP , NH  ,FT ,FD , FP , FH 
GO  TO  (5, 6, 7, 8) ,IT 

5 TTI =T 
RETURN 

6 TTI =T-273 . 15D0 
RETURN 

7 TTI=T*1.8D0 
RETURN 

8 TTI=T*1.6D0-45S.67D0 
RETURN 

END 


APPENDIX  B 


This  appendix  contains  sample  calculations  to  help  the  user  test 
his  program.  Also  listed  are  the  main  programs  used  to  print  out 
numerical  values. 


C THIS  IS  A MAIN  PROGRAM  FOR  THE  CALCULATION  OF  TABLES  OF  PROPERTIES 
C USING  A CHOICE  OF  UNITS,  AND  AT  A CHOICE  OF  CONSTANT  TEMPERATURE, 

C PRESSURE  OR  DENSITY.  THE  USER  IS  QUERIED  AS  TO  HIS  CHOICES  OF 
C UNITS  AND  VARIABLES  AT  EXECUTION  TIME. 

IMPLICIT  DOUILE  PRECISION ( A-H , O-Z) 

DOUBLE  PRECISION  NT  ,ND , NP , NH 

COMMON  /UNITS/  IT , ID , IP  , IH , NT , ND , NP  , NH ,FT ,FD ,FP ,FH 
COMMON  /QQQQ/  Q0,Q5 

COMMON  /FCTS/  AD , GD  ,SD , UD ,HD , C VD , CPD ,DPDT 
COMMON  /ACONST/  WM, GASCON, TC,AA,Z,DZ,Y  .UREF  ,SREF 
COMMON  /N CONST / G (40  ) , 1 1 (40 ) , J J ( 40 ) ,NC 
DATA  NS1  ,NS2/2H  m,2Hf t/ 

CALL  UNIT 
NS=NS1 

IF ( ID  .EQ  .4  > NS  = NS2 
15  WRITE (6, 11 ) 

READ ( 5 ,* ,£ND=9  ) IOPT.XISO 
IF ( IOPT . EQ  .0  ) GO  TO  9 
GO  TO  (101,201,301)  , IOPT 

101  WRITE(6, 102  ) 

102  FORMAT ( ' ENTER  OTHER  INDEPENDENT  VARIABLE  (2  FOR  P,  3 FOR  DENS'/ 

1'  FOLLOWED  BY  INITIAL,  FINAL  AND  INCREMENTAL  VALUES  OF  THIS  VAR.') 

READ(5,*,END=9)  JOPT, Y1 , Y2 , YI 
IF(JOPT-i)  15,101,103 

103  TT=XI SO 
T=TTT  (TT  ) 

WRITE (6, 444)  NT ,NP , ND , NP , NT ,NH  ,NT , NH ,NS 

IF( JOPT.EC  .2)  DGSS=Yl/FP/T/.4 

IZ=0 

CALL  3B( T ) 

P3S=20000 . 

DVV=0  . 

IF(T.LT.TC)  CALL  PCORR( T,PSS, DLL ,DVV ) 

DGSS=DVV 

IF( DGSS. EQ  .0  . ) DGSS=1.11-.0004*T 
PSAT=PSS*FP 

I F ( JOPT . EQ  .2  .AND.  Yl.GT.PSAT)  IZ=3 

IF(Yl.GT.PSAT)  DGSS=DLL 

IF ( JOPT. GE .3 ) I Z=3 

FI N=Y 1-Y I 

DIN=PIN 

PINC=YI/FP 

DINC=YI*FD 

22  IF( JOPT.EQ  .2)  PIN=PIN+YI 
II( JOPT. EQ. 3)  DIN-DI N+YI 
I F( JOPT . EQ  .2  .AND.  PIN.GT.Y2)  GO  TO  101 
IF( JOPT. EQ  .3  .AND.  DIN.GT.Y2)  GO  TO  101 
I F (JOPT . EQ .2 ) PRES  =P IN/FP 
I F ( JOPT. EC. 3)  D-DIN/FD 
24  CONTINUE 


23 


26 

27 


IF( JOPT.EQ .2 
IKJOPT.EQ.1 
TSAYE=TT-YI 
IF ( JOPT . EQ . 1 
PSAVE=PI N -YI 
IF( JOPT.EQ.  2 
IF( JOPT. EQ .1 
IZ=IZ+1 


.OR.  (JOPT.EQ. 2 .AND.  PIN  .LT  .PSAT ) ) GO  TO  26 
.AND.  T.LT.TS)  GO  TO  26 

.AND.  I0PT.EQ.2  .AND.  IZ.LE.2)  TT=TT I (TS ) 

.AND.  PRES .GT.PSAT/FP  .AND.  IZ.GE.2)  GO  TO  26 
.AND.  T.GT.TS  .AND.  IZ.GE.2)  GO  TO  26 


IF( JOPT.EQ. 2)  PRES  =PSAT /FP 
IF ( JOPT.EQ.  1)  T=TS 

IF ( I Z . EQ  . 1 .AND.  JOPT.EQ.l)  DG3S=DLL 
IF(IZ.EQ.l  .AND.  JOPT.EQ. 2)  DG3S=DYY 
IF( IZ  .EQ .2  .AND.  JOPT.EQ. 2)  DGSS=DLL 
IF ( I Z .EQ . 2 .AND.  JOPT.EQ.l)  DGSS=DYY 
C AL  L E B ( T ) 

IF ( IOPT . NS  .3  .AND.  J0PT.NE.3)  CALL  DFI ND (D ,PRES , DGSS , T , DQ ) 
CALL  QQ( T ,D  ) 

CALL  BA5E(D ,T ) 

RT=GASCON*T 


PDUM  = D*RT*Z  + Q0 

IF ( JOPT . EQ  .3  .OR.  I0PT.EQ.3)  PRES=PDUM 
I F ( IOPT . EQ  .3  .OR.  JOPT.EQ. 3)  DQ=RT*( Z+Y*D Z ) +Q5 
DGSS=D  + PINC/DQ 
CALL  THERM(D ,T ) 

U = UD*T*GASCON*FH 


C=DSQRT( DABS ( CPD*DQ*1 . D3/CYD ) ) 

IF( ID .EQ .4 ) C=C*3 . 280833 

H = HD^T^GASCON^FH 

S = SD*GASCON*FH*FT 

dpdti=dpdt:J:f  p*f  t 

d'odd=d%*f  d 

COMP  = 1 .E3/D/DQ/FP 

DDDTL=1.D3*DPDT/D/DQ 

CP=CPD*GASCON*FH*FT 


CV=CVD*GASCON*FH*FT 

VL=fd/D 


DOUT=l./YL 
’ POUT=PRES’!‘FP 

WRITE  (6, 21)  TT,POUT,DOUT,DPDTX tDFDD,CY  ,CPfS  ,H,U,C 
21  FORMAT  (F9 .3  ,F12 .5  , F12. 7 t FH  • 5,  FI  1 .3 , 5F12  .4  , Fll  .3) 

444  FORMAT ( 5X , 1ST , 12X , 1HP , 10X , 1HD ,8X , 5HdP/dT , 5X , 5HdP/dD , 10X , 2HCr 

l,10Xf 2HCp,llXf 1HS ,9X , 1HH , 11X , 1HU ,7X  f 7HYel  Snd/3X,4Hdeg  ,A1,7X,A8 

2 , 4X ♦ A6 , 5X  ,A3 ♦ 1H/ , A1 ,4X ♦ 7H  ,15X,4H ,A6,A1,4H 

3 , 13X, 4H ♦ A6 ,4H  - - ,8X ,A2 , 4H/sec/ ) 

11  FORMAT ( ' DO  YOU  WISH  TO  CALCULATE  AN  ISOTHERM  (ENTER  1),  AN  ISOBA 
1R  (ENTER  2)  OR  AN  ISOCHORE?  (ENTER  3)V'  FOLLOWED  BY  YALUS  OF  ISO. 
2 (ENTER  0 TO  DISCONTINUE)') 

IF(IZ.EQ.l)  WRI TE ( 6 , 12 ) 

12  FORMATd  


) 


IF(IZ.EQ.l)  GO  TO  22 

IF(IZ.EQ.2  .AND.  J0PT.EQ.2)  PIN=FSAVE 
IF(IZ.EQ.2  .AND.  JOPT.EQ.l)  TT=TSAVE 
IE ( IZ . EQ . 2 ) IZ  = 3 
GO  TO  (22,204,304 ) ,IOPT 

201  JOPT=l 

PRES  = XISO/F? 

202  WR I TE ( 6 203) 

203  FORMAT ( ' ENTER  FIRST  LAST  AND  INCREMENT  OF  T ' ) 
READ( 5 ,END=S ) T1,T2,YI 

IF(T1.EQ.0.)  GO  TO  15 

WRITE (6, 444)  NT  ,NP  ,ND , NP  ,NT  ,NH , NT , NH  ,NS 

TT=T1-YI 

T=TTT (TT ) 

CALL  TCORR( TS ,PRES  fDLL ,DVV ) 

D = DLL 

IF( (T+YI*FT) .GT .TS ) D=DVV 
I Z=3 

I F ( (T+YI *FT ) .LT .TS )IZ=0 

204  TT=TT+YI 
T=TTT (TT ) 

IF(TT.GT.T2)  GO  TO  202 
CALL  BB(T) 

DGS  S=D 
GO  TO  24 
201  JOPT=l 
D=XI SO 

302  WRITE (6, 203 ) 

READ ( 5 , END=9 ) T1,T2,YI 
IF ( T1 .LE  .0  . ) GO  TO  15 

WRITE (6, 444)  NT ,NP  ,ND , NP , NT  , NH  ,NT  , NH  ,NS 
TT=T1-YI 
I Z=3 

T=TTT( TT) 

204  TT=TT+YI 
T=TTT ( TT ) 

IF( TT . GT . T2 ) GO  TO  302 
CALL  BB(T) 

GO  TO  27 
9 STOP 

END 


o o o o a o 


THE  FOLLOWING  IS  A MAIN  PROGRAM  FOR  CALCULATING  THE  PROPERTIES 

ALONG  THE  SATURATION  CURVE  UP  TO  646.3  I.  THE  INPUTS  ARE:  INITIAL  T, 
FINAL  T AND  THE  T INCREMENT  (ALL  IN  K).  THE  VAPOR  PRESSURE  IS  FIRST 
CALCULATED  USING  "PS"  AND  ADJUSTED  FROM  THERE  SO  THAT 
GL-GV=0.  ALSO  IN  THE  INPUT  ARE  AN  INITIAL  GUESS  FOR  THE  LIQUID  AND 
VAPOR  DENSITIES.  IF  NOT  GIVEN,  THEY  WILL  BE  INTERNALLY  CHOSEN. 
IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /ACONST/  WM .GASCON, TZ ,AA, ZB ,DZB ,YB 
COMMON  /FCTS/  AD , GD ,SD ,UD ,HD , CVD  ,CPD ,DPDT 

3 READ ( 5 , 4, END =99  ) TI  ,TF,TD,DL  ,DG 

4 FORMAT () 

WRITE ( 6,22) 

T=TI-TD 

6 T=T+TD 

IF(T.GT.TF)  GO  TO  3 
IF (T  ,GT  .646.3D0)  GO  TO  3 
I JK=0 

CALL  BB (T ) 

P=PS ( T)  +1 JK*PT 

C THE  FOLLOWING  TWO  LINES  CALCULATE  AN  INITIAL  GUESS  FOR  THE  DENSITIES 
C IF  GUESSES  NOT  SUPPLIED  IN  THE  INPUT. 

IF (DL . LE.0 . ) DL=1 . 11- . 0004*T 
IF(DG .LE.0 . ) DG=P/GASCON/T 

5 P=P+IJK*PT 
RT=GASCON*T 

CALL  DFIND (RL ,P , DL , T, DPDL) 

DL=RL 

CALL  THERM (RL ,T ) 

SL=SD*GASCON 
HL  = HD^GASCON^T 
VL=1./DL 
GL=GD 

CALL  DFIND (RG ,P ,DG , T, DPDG) 

DG=RG 

IF (RG .LT .5 .D-7 ) RG=5.D-7 
CALL  THERM ( RG  ,T  ) 

SG=SD*GASCON 
EG  = HD*GASCON*T 
VG=1./DG 
GV=GD 

HEAT=HG-HL 

DELG=GL-GV 

I F ( DABS (DELG ) . LT .2 .D-6 ) GOTO  15 
PT  = DELG^GASCON^T/ ( VG-VL ) 

IF(T  .GT  .640.D0 ) PT=.1*PT 
I JE=I JK+1 

IF(IJK.EQ.l  .OR.  (T.GT.64.D1  .AND.  IJI.LE.10))  GO  TO  5 


15 

HEAT=HG-HL 

WRITE (6,21)  T , P ,RL 

, RG ,HL ,HG  fHEAT ,SL  ,SG,VL,VG 

21 

FORMAT ( F9.3,F13 .6, 

F9.6  ,F9.7 ,3FS  .2  ,2F9 .4,F9 .5,F11 .3) 

22 

FORMAT ( ' T,K 

P.MPA  DL  ,G/CC  DV.G/CC  HL,J/G 

1 LAT  HT  SL.J/GK 

3V,J/G5  VL  VG ' /) 

GO  TO  6 

99 

STOP 

END 

Ef> 


C THIS  IS  A SAMPLE  MAIN  PROGRAM  WHICH  WILL  SERVE  AS  AN  EXAMPLE  FOR 
C THE  USE  OF  THE  SUBROUTINES  AND  FUNCTIONS  GIVEN  ABOVE,  AND  WHICH 
C WILL  PRINT  OUT  VALUES  OF  VARIOUS  PROPERTIES  CALCULATED  FOR  A 
C GIVEN  INPUT  POINT  TO  A LARGE  NUMBER  OF  SIGNIFICANT  FIGURES, 

C SUITABLE  FOR  USE  AS  A CHECK  ON  THE  OPERATION  OF  THE  PROGRAM. 

C THE  USER  IS  QUERIED  AS  TO  THE  UNITS  DESIRED. 

IMPLICIT  DOUBLE  PRECISION (A-H ,0-Z ) 

DOUBLE  PRECISION  NT ,ND , NP , NH 

COMMON  /UNITS/  IT , ID , IP , IH , NT , ND , NP , NH ,FT  ,FD , FP , FH 
COMMON  /QQQQ/  Q0,Q 5 

COMMON  /FCTS/  AD  ,GD  ,SD ,UD , HD , C VD  , CPD  ,DPDT ,DVDT , C JTT , CJTH 
COMMON  /ACONST/  WM , GASCON , TZ ,AA, Z ,DZ ,Y ,UREF ,SREF 
COMMON  /NCONST/  G ( 40  ) , 1 1 (40  ) , J J ( 40  ) , NC 
DATA  NSS1,NSS2/2H  m,2Hft/ 

CALL  UNIT 
N S=NS  S 1 

IF( ID . EC .4 ) N S=NSS 2 
15  WRITE(6, 11) 

100  READ( 5,* ,END=9)  IOPT,X,TT 
T=TTT( TT ) 

RT=GAS  CO  N*T 
CALL  EB( T ) 

IF( IOPT.LE .0 ) GOTO  9 

GOTO  (101, 102, 100, 100, 100), IOPT 

101  DD=X 
D=DD*FD 
CALL  CO ( T ,D ) 

CALL  BASE (D ,T ) 

PRES  = FP* ( RT*D*Z  + 00) 

DQ=RT* ( Z+Y*DZ ) +05 
GOTO  111 

102  PRES=X 
P=PRES/FP 
DGS  S=P/T / .4 
PSAT=20000. 

IF(T.LT.TZ)  CALL  PCORR (T  , PS AT, DLL ,DGSS ) 

I F ( P . GT . PSAT ) DGSS=DLL 

IF( P. GE. PSAT  .AND.  T.LT.523.15)  DGSS=1 . 11-. 0004*T 
CALL  DFI ND(D,P,DGSS,T,DQ) 

DD=D/FD 

111  CALL  THERM ( D , T ) 

U = UD*T*GAS  CON*FH 

C=DSORT ( DABS ( CPD*D0*1 . D3/CVD ) ) 

IF ( ID  . EO .4 ) C=C*3. 280833 
H = HD*T*GASCON*FH 
S = SD*GASCON*FH*FT 
CP=C?D*GASCON*FH*FT 
C V=C VD*G ASCON*FH*FT 
VL=1./D 

DPDD  = D0’!tFD*FP 
DPDT1=DPDT*FP*FT 

WRI TE ( 6 , 20 ) TT  ,NT  ,PRES , NP ,DD, ND  j 


£7 


WRITE(e,21)  DPDTl,DPDD,CT,NH,NT,CP,S,H,NH,n,C,NS,CJTT,CJTH,DVDT 

20  FORMAT( ' T = ,F12.4,  deg  ,A1 ,5X»  P = »F13 .6, IX » A6, 5X * D = 
1,F14.10,1X,A8) 

21  FORMAT ( ' DP/DT  = ' ,F16 .9 ,6X , 'DP/DD  =',F16.5/ 

2 ' C V =',F12.6,1X,A6,A1,5X,'CP  = ' ,F12  .6  ,6X , 'S  =',F12.6/ 

3 ' H =',F14.6,1X,A6,5X, 'U  = ' ,F14 .6 ,6X , 'VEL  SND  = ' , F14. 6 ,A2 , '/sec '/ 

4 ' JT( T)  =' ,F11 .5 ,5X, 'JT(H)  = ' , F 11 .5 ,5X , 'DV/DT  =',F12.6/) 

11  FORMAT ( ' INTER  OPTION,  X,  AND  T,  WHERE  FOR  OPTION=l,  X=DENSITY'/ 

1'  AND  FOR  0PTI0N=2 , X=PRESSURE  (ENTER  0 TO  QUIT)')  . 

GO  TO  100 
9 STOP 

END 


✓jXQT  POINT 

ENTER  UNITS  CHOSEN  FOR  TEMPERAT 

CHOOSE  FROM  1=DEG  X,  2=DEG  C,  3=D£G  R,  4=DEG  F 

>1 

ENTER  UNITS  CHOSEN  FOR  DENSITY 
CHOOSI  FROM  1=KG/M3 , 2=G/CM3,  3=M0L/L,  4=LB/FT3 
>2 

ENTER  UNITS  CHOSEN  FOR  PRESSURE  • 

CHOOSE  FROM  1=MPA  , 2=BAR , 3=ATM,  4=PSIA,  5=XG/CM2 

>2 

ENTER  UNITS  CHOSEN  FOR  ENERGY 

CHOOSE  FROM  1=X J/XG , 2=J/G,  3=J/M0L,  4=CALORIES /G , 5=CAL0RI ES/MOL , 6=BTU/L3 

>2 

ENTER  OPTION,  X,  AND  T,  WHERE  FOR  OPTION=l,  X=DENSITY 
AND  FOR  OPTI 0N=2 , X=PRESSURE  (ENTER  2 TO  QUIT) 

>1,  .9,  873.15 

T = 673.1500  deg  X P = 7110.605028  Bar  D = .9000000000  g/cm3 

DP/DT  = 14.491600834  DP/DD  = 28719.49752 

C7  = 2.827220  J/g  X CP  = 3.615462  5 = 4.064690 

H = 2779.151751  J/g  U = 1989.062303  VEL  SND  = 1916.419293  m/sec 

JT( T ) = .56718  JT(H)  = -.15688  DV/DT  = .000623  - 


>2  , 225.  , 648.15 
T = 648.1500  deg  K 

DP/DT  = 3.447886683 

CV  = 3.743787  J/g  X 

H = 1965.692198  J/g 


P = 225.000000  Bar  D = .4103745556  g/cm3 

DP/DD  = 63.95378 

CP  = 75.284775  S = 4.221788 

U = 1910.664237  VEL  SND  = 358.617190  m/sec 


JT ( T ) = -205.05536 


JT(H)  = 


2.72373 


DV/DT  = 


.320130 


>2,  .00617,  273.16 

T = 273.1600  deg  K 

DP/DT  = -1.576872063 

CV  = 4.225225  J/g  K 

H = .000617  J/g 


JT  ( T ) = 


P = .006170  Bar  D = .9997782169  g/cm3 

DP/DD  = 19608  .40085 

CP  = 4.228690  S = .000000 

U = -.000000  VEL  SND  = 1400.674132  n/sec 


1 .02220 


JT (H ) = 


-.24173 


DV/DT  = 


- .000080 


> 


T 

P 

D 

V 

deg  C 

MPa 

g/cm3 

cm3/g 

50.000 

.0100 

. 000067 

14869.238335 

50.000 

.0123 

.000083 

12036.663719 

50.000 

.0123 

.987991 

1.012155 

50.000 

.1000 

.988030 

1 .012116 

50.000 

1 .0000 

.988422 

1.011713 

50 .000 

10 .0000 

.992305 

1.007754 

50.000 

20 . 0000 

.996526 

1.003486 

50 .000 

30.0000 

1 . 000656 

.999344 

50.000 

40 .0000 

1.004700 

.995322 

50.000 

50 .0000 

1.008663 

.991411 

50.000 

60  . 0000 

1.012550 

.987605 

50.000 

70  .0000 

1.016365 

.963899 

50 .000 

80.0000 

1.020109 

.980287 

50.000 

90 . 0000 

1.023788 

.976765 

50.000 

100  .0000 

1.027403 

.973328 

50.000 

200 . 0000 

1.060447 

.942999 

50.000 

300.0000 

1.088794 

.918446 

50 . 000 

400  . 0000 

1.113570 

.898013 

50 .000 

500.0000 

1.135810 

.880429 

50 .000 

600 . 0000 

1.156346 

.864793 

50  .000 

700  .0000 

1.175655 

.850590 

50.000 

800 .0000 

1.193766 

.837685 

50 .000 

900  .0000 

1.210356 

.826203 

50 .000 

1000  .0000 

1.225100 

.816260 

T 

P 

D 

V 

deg  C 

MPa 

ta  / c m3 

cm3/g 

250.000 

.0100 

. 000041 

24136.161956 

250.000 

.1000 

. 000416 

2406.053164 

250 .000 

1 .0000 

. 004298 

232.644926 

250.000 

3.9736 

.019956 

£0 .110581 

250 .000 

3 .9736 

.799072 

1.251452 

250 .000 

10  .0000 

.805899 

1.240850 

250 . 000 

20  .0000 

.816284 

1.225064 

250 .000 

30  . 0000 

.825733 

1.211046 

250 .000 

40  .0000 

.834438 

1.198411 

250 .000 

50  .0000 

.842537 

1.186892 

250  .000 

60  .0000 

.350128 

1.176293 

250 .000 

70  .0000 

.857288 

1.166469 

250  .000 

80  .0000 

.864075 

1.157307 

250.000 

90  .0000 

.870537 

1.148717 

250.000 

100  .0000 

.876711 

1.140627 

252 . 000 

200  .0000 

.927623 

1 .078024 

250 .0  00 

300  .0000 

.966627 

1.034311 

250.002 

400  .0000 

.999498 

1 .000502 

250.000 

500  .0000 

1.027884 

.972872 

250.000 

600 .0000 

1 .053194 

.949493 

250.000 

700  .0000 

1.076161 

.929229 

250 . 000 

800  .0000 

1.097271 

.911352 

250.000 

900 .0000 

1.116861 

.895367 

250 . 020 

1000  .0000 

1 . 135179 

.880918 

T 

P 

D 

V . 

deg  C 

MPa 

g /cm3 

cm3/g 

375.000 

.0100 

.000033 

29909.055417 

375.000 

.1000 

.000335 

2986.855135 

375.000 

1 .0000 

.003395 

294.567937 

375.000 

10  . 0000 

.040763 

24.532263 

375.000 

20  . 0000 

.130420 

7.667535 

375.000 

22  .3292 

.320009 

3.124909 

375.000 

30 . 0000 

.558254 

1.791298 

375.000 

40 . 0000 

.609556 

1.640537 

375.000 

50.0000 

.641323 

1.559276 

375.000 

50  .0000 

.665164 

1.503388 

375.000 

70  .0000 

.664567 

1.460777 

375.000 

80 . 0000 

.701090 

1.426350 

375.000 

90.0000 

.715576 

1 .397475 

375.000 

100  . 0000 

.728535 

1.372617 

375 . 000 

200 .0000 

.816920 

1.224110 

375.000 

300 .0000 

.873641 

1.144634 

375.000 

400  .0000 

.917244 

1.090222 

375.000 

500  .0000 

.953396 

1 .048882 

375.000 

600 . 0000 

.984657 

1.0155e2 

375.000 

700 .0000  e 

1.012420 

.987732 

375.000 

800 . 0000 

1.037537 

.963821 

375.000 

900  .0000 

1.060570 

.942889 

375.000 

1000  .0000 

1.081910 

.924291 

T 

P 

D 

V 

deg  C 

MPa 

g /cm3 

cm3/g 

500.000 

.0100 

.000028 

35679.867768 

500.000 

.1000 

.00 028 e 

3565.549685 

500 . 000 

1.0000 

.002824 

354.099874 

500 .000 

10.0000 

. 030503 

32.783626 

500 .000 

20 .0000 

. 067711 

14.768678 

500  .000 

30 .0000 

.115259 

8.676122 

500.000 

40.0000 

. 177974 

5.618801 

500 . 000 

50  .0000 

.256947 

3.891860 

500 .000 

60 . 0000 

.338443 

2.954707 

50 0 . 000 

70  .0000 

.405762 

2.464496 

500 .000 

80  .0000 

.456993 

2.188216 

500 . 000 

90  .0000 

.496532 

2.013967 

500.000 

100  .0000 

.528211 

1 .893181 

500  .000 

200  .0000 

.691177 

1.446808 

500 .000 

300  . 0000 

.772885 

1.293853 

500 . 000 

400  .0000 

.830456 

1 .204158 

500 .000 

500  .0000 

.876029 

1.141515 

500  .000 

600  .0000 

.914303 

1.093730 

500 .000 

700  .0000 

.947615 

1 .055281 

500 .000 

800  . 000 0 

.977308 

1.023219 

500 .000 

900  .0000 

1.004227 

.995791 

50  0 . 000 

1000  .0000 

1.028941 

.971873 

■\ 


B 1° 


T 

P 

D 

V 

deg  C 

MPa 

S /cm3 

cm3/g 

750 .000 

.0100 

.000021 

47219.519942 

750 . 000 

.1000 

.000212 

4720 .958893 

750.000 

1 .0000 

.002123 

471.103828 

750.000 

10  . 0000 

.021679 

46.127641 

750.000 

20  .0000 

. 044390 

22.527645 

750.000 

30  .0000 

.068163 

14.670731 

750 .000 

40  .0000 

.092985 

10.754408 

750 .000 

50  .0000 

.118770 

8.419521 

750 . 000 

60  .0000 

.145340 

6.880435 

750.000 

70  .0000 

. 172425 

5.799609 

750 . 000 

80.0000 

. 199697 

5.007588 

750.000 

90  .0000 

.226804 

4.409089 

750 .000 

100  .0000 

.253421 

3.946003 

750.000 

200  .0000 

.462110 

2.163988 

750 .000 

300 .0000 

.584047 

1.712190 

750 . 000 

400  .0000 

.665974 

1.501561 

750.000 

500  .0000 

.728216 

1.373220 

750.000 

600  .0000 

.779010 

1.283680 

750 . 000 

700  .0000 

.822326 

1.216062 

750 .000 

800  .0000 

.860356 

1.162309 

750 .000 

900  . 0000 

.894435 

1.118025 

750.000 

1000  .0000 

. .925435 

1.080573 

T 

P 

D 

V 

deg  C 

MPa 

g 'cm3- 

cm3/g 

1000.000 

.0100 

.000017 

58758.280969 

1000.000 

.1000 

.000170 

5875.475023 

1000 . 000 

1 .0000 

.001703 

587.198410 

1000.000 

10 . 0000 

.017122 

58.404407 

1000 . 000 

20 . 0000 

.034420 

29.052786 

1000 .000 

30 . 0000 

.051866 

19  .280576 

1000 . 000 

40 . 0000 

.069439 

14.401171 

1000. 000 

50  .0000 

.087117 

11 .478860 

1000.000 

60  .0000 

.104867 

9.535858 

1000.000 

70.0000 

. 122648 

6.153425 

1000.000 

£0 . 0000 

. 140405 

7.122229 

1000.000 

90  .0000 

. 156081 

6 .325872 

1000.000 

100  .0000 

. 175613 

5.694349 

1000.000 

200  .0000 

.333285 

3.000438 

1000.000 

300 . 0000 

.451565 

2.214522 

1000.000 

400  .0000 

.539679 

1.652267 

1000.000 

500 . 0000 

.609226 

1.641425 

1000 .000 

600  .0000 

.666407 

1.500564 

1000 . 000 

700 . 0000 

.715336 

1 .397941 

1000.000 

800 . 0000 

.756357 

1.318641 

1000.000 

900  .0000 

.796939 

1.254601 

1000 .000 

1000 .0000 

.832065 

1 .201830 
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